# All functions needed to run network extension 

## Start out by defining parameter structure 
@with_kw struct define_param_industry
    connect_theta::Float64 = 0.2 
    connect_cost_level::Float64 = 1.014782e-08 
    connect_cost_elasticity::Float64 = 1.04
    connect_fixed_cost::Float64 = 0.0
    connect_rho::Float64 = -0.999 
    alpha_tilde::Float64 = 0.15
    beta_tilde::Float64 = 0.2
    gamma_tilde::Float64 = 0.5
    eta_tilde::Float64 = 0.85
    elasticity_ratio::Float64 = 0.85/(1-0.85)
    sigma::Float64 = 1/(1-0.85)
    alpha_eps::Float64 = 7.016288e-02
    beta_eps::Float64 = -1.311290e-02
    variance_eps_z::Float64 = 5.353032e-07
    proba_c::Float64 = baseline_parameters_df.proba_c[1]
    share_industry::Float64 = 1/6
end

# want to be able to index the structure
Base.getindex(x::define_param_industry, i) = getproperty(x, fieldnames(typeof(x))[i])

# Read z_star_grids (should be arrays of the form: [z,j])
z_star_grid_network = DataFrame(CSV.File(joinpath(ddat,"industry_z_star_grid.csv"); header=1))[:,:]
z_star_distrib_network = DataFrame(CSV.File(joinpath(ddat,"industry_z_star_distribution.csv"); header=1))[:,:]

# save number of sectors 
const n_sectors = 6::Int64

# get parameters for chemicals 
chemicals_parameters = DataFrame(CSV.File(joinpath(ddat,"industry_parameters_chemicals.csv"); header=1))
chemicals_parameters = define_param_industry(
    connect_theta = chemicals_parameters.theta[1],
    connect_cost_level = chemicals_parameters.cost_level[1],
    connect_cost_elasticity = chemicals_parameters.cost_elasticity[1],
    connect_fixed_cost = 0.0, 
    connect_rho = chemicals_parameters.rho[1],
    alpha_tilde = chemicals_parameters.alpha_tilde[1],
    beta_tilde = chemicals_parameters.beta_tilde[1],
    gamma_tilde = chemicals_parameters.gamma_tilde[1],
    eta_tilde = chemicals_parameters.eta_tilde[1],
    elasticity_ratio = chemicals_parameters.eta_tilde[1]/(1-chemicals_parameters.eta_tilde[1]),
    sigma = 1/(1-chemicals_parameters.eta_tilde[1]),
    alpha_eps = chemicals_parameters.alpha_eps[1],
    beta_eps = chemicals_parameters.beta_eps[1],
    variance_eps_z = chemicals_parameters.variance_eps_z[1],
    proba_c = chemicals_parameters.proba_c[1], 
    share_industry = chemicals_parameters.share_industry[1]) 

# get parameters for machinery
machinery_parameters = DataFrame(CSV.File(joinpath(ddat,"industry_parameters_machinery.csv"); header=1))
machinery_parameters = define_param_industry(
    connect_theta = machinery_parameters.theta[1],
    connect_cost_level = machinery_parameters.cost_level[1],
    connect_cost_elasticity = machinery_parameters.cost_elasticity[1],
    connect_fixed_cost = 0.0, 
    connect_rho = machinery_parameters.rho[1],
    alpha_tilde = machinery_parameters.alpha_tilde[1],
    beta_tilde = machinery_parameters.beta_tilde[1],
    gamma_tilde = machinery_parameters.gamma_tilde[1],
    eta_tilde = machinery_parameters.eta_tilde[1],
    elasticity_ratio = machinery_parameters.eta_tilde[1]/(1-machinery_parameters.eta_tilde[1]),
    sigma = 1/(1-machinery_parameters.eta_tilde[1]),
    alpha_eps = machinery_parameters.alpha_eps[1],
    beta_eps = machinery_parameters.beta_eps[1],
    variance_eps_z = machinery_parameters.variance_eps_z[1],
    proba_c = machinery_parameters.proba_c[1], 
    share_industry = machinery_parameters.share_industry[1]) 

# get parameters for food 
food_parameters = DataFrame(CSV.File(joinpath(ddat,"industry_parameters_food.csv"); header=1))
food_parameters = define_param_industry(
    connect_theta = food_parameters.theta[1],
    connect_cost_level = food_parameters.cost_level[1],
    connect_cost_elasticity = food_parameters.cost_elasticity[1],
    connect_fixed_cost = 0.0, 
    connect_rho = food_parameters.rho[1],
    alpha_tilde = food_parameters.alpha_tilde[1],
    beta_tilde = food_parameters.beta_tilde[1],
    gamma_tilde = food_parameters.gamma_tilde[1],
    eta_tilde = food_parameters.eta_tilde[1],
    elasticity_ratio = food_parameters.eta_tilde[1]/(1-food_parameters.eta_tilde[1]),
    sigma = 1/(1-food_parameters.eta_tilde[1]),
    alpha_eps = food_parameters.alpha_eps[1],
    beta_eps = food_parameters.beta_eps[1],
    variance_eps_z = food_parameters.variance_eps_z[1],
    proba_c = food_parameters.proba_c[1], 
    share_industry = food_parameters.share_industry[1]) 

# get parameters for minerals  
minerals_parameters = DataFrame(CSV.File(joinpath(ddat,"industry_parameters_minerals.csv"); header=1))
minerals_parameters = define_param_industry(
    connect_theta = minerals_parameters.theta[1],
    connect_cost_level = minerals_parameters.cost_level[1],
    connect_cost_elasticity = minerals_parameters.cost_elasticity[1],
    connect_fixed_cost = 0.0, 
    connect_rho = minerals_parameters.rho[1],
    alpha_tilde = minerals_parameters.alpha_tilde[1],
    beta_tilde = minerals_parameters.beta_tilde[1],
    gamma_tilde = minerals_parameters.gamma_tilde[1],
    eta_tilde = minerals_parameters.eta_tilde[1],
    elasticity_ratio = minerals_parameters.eta_tilde[1]/(1-minerals_parameters.eta_tilde[1]),
    sigma = 1/(1-minerals_parameters.eta_tilde[1]),
    alpha_eps = minerals_parameters.alpha_eps[1],
    beta_eps = minerals_parameters.beta_eps[1],
    variance_eps_z = minerals_parameters.variance_eps_z[1],
    proba_c = minerals_parameters.proba_c[1], 
    share_industry = minerals_parameters.share_industry[1]) 

# get parameters for wood 
wood_parameters = DataFrame(CSV.File(joinpath(ddat,"industry_parameters_wood.csv"); header=1))
wood_parameters = define_param_industry(
    connect_theta = wood_parameters.theta[1],
    connect_cost_level = wood_parameters.cost_level[1],
    connect_cost_elasticity = wood_parameters.cost_elasticity[1],
    connect_fixed_cost = 0.0, 
    connect_rho = wood_parameters.rho[1],
    alpha_tilde = wood_parameters.alpha_tilde[1],
    beta_tilde = wood_parameters.beta_tilde[1],
    gamma_tilde = wood_parameters.gamma_tilde[1],
    eta_tilde = wood_parameters.eta_tilde[1],
    elasticity_ratio = wood_parameters.eta_tilde[1]/(1-wood_parameters.eta_tilde[1]),
    sigma = 1/(1-wood_parameters.eta_tilde[1]),
    alpha_eps = wood_parameters.alpha_eps[1],
    beta_eps = wood_parameters.beta_eps[1],
    variance_eps_z = wood_parameters.variance_eps_z[1],
    proba_c = wood_parameters.proba_c[1], 
    share_industry = wood_parameters.share_industry[1]) 

# get parameters for textiles 
textiles_parameters = DataFrame(CSV.File(joinpath(ddat,"industry_parameters_textiles.csv"); header=1))
textiles_parameters = define_param_industry(
    connect_theta = textiles_parameters.theta[1],
    connect_cost_level = textiles_parameters.cost_level[1],
    connect_cost_elasticity = textiles_parameters.cost_elasticity[1],
    connect_fixed_cost = 0.0, 
    connect_rho = textiles_parameters.rho[1],
    alpha_tilde = textiles_parameters.alpha_tilde[1],
    beta_tilde = textiles_parameters.beta_tilde[1],
    gamma_tilde = textiles_parameters.gamma_tilde[1],
    eta_tilde = textiles_parameters.eta_tilde[1],
    elasticity_ratio = textiles_parameters.eta_tilde[1]/(1-textiles_parameters.eta_tilde[1]),
    sigma = 1/(1-textiles_parameters.eta_tilde[1]),
    alpha_eps = textiles_parameters.alpha_eps[1],
    beta_eps = textiles_parameters.beta_eps[1],
    variance_eps_z = textiles_parameters.variance_eps_z[1],
    proba_c = textiles_parameters.proba_c[1], 
    share_industry = textiles_parameters.share_industry[1]) 

## Create network_baseline_objects
network_baseline_objects = (
    z_star_grid = Matrix(z_star_grid_network), 
    z_star_distrib = Matrix(z_star_distrib_network), 
    param = [
        chemicals_parameters, 
        machinery_parameters,
        food_parameters,
        minerals_parameters,
        wood_parameters,
        textiles_parameters]
)

### Get further network elasticities 

# final use shares 
const network_nu = [
    0.0597145, # chemicals 
    0.1852055, # machinery
    0.504857, # food 
    0.004605, # minerals 
    0.023186, # wood 
    0.222432 # textiles 
     ]

# Get elasticities 
industry_io_elasticities = (Matrix(DataFrame(CSV.File(joinpath(ddat,"industry_io_elasticities.csv"); header=1))))'

## Make sure that industry IO elasticities are correct: 
industry_io_elasticities = industry_io_elasticities ./ sum(industry_io_elasticities,dims=2)

# estimate subsidy for one industry over fully specified grid 
function compute_subsidy_industry(; z_star_grid, epsilon_grid, param)

    #### Step 1: Preparation ####

    # get parameters 
    theta = param.connect_theta 
    cost_level = param.connect_cost_level
    cost_elast = param.connect_cost_elasticity 
    fixed_cost = param.connect_fixed_cost 

    # baseline x_star (enforcing w = 1.0 & leaving out P_j_m (also left out in z_star, so cancels out)) 
    x_star = ( (((1-vat_tax)*param.alpha_tilde/(rental_rate - 1.0))^param.alpha_tilde)*(((1-vat_tax)*param.beta_tilde)^param.beta_tilde)*((param.gamma_tilde)^param.gamma_tilde) )

    # get z_tilde_grid 
    z_tilde_grid = (z_star_grid .* x_star).^(1/(1-param.eta_tilde))
  
    # specify FOC (this relies on having m_R paying P, not P_j_m)
    function rent_seeking_FOC(m_R, z_tilde_value, epsilon)
        return z_tilde_value .* ( (1.0 .+ (epsilon .*((m_R).^(theta)) .- cost_level .* (m_R.^cost_elast))).^param.elasticity_ratio ) .* ( (theta .* epsilon .* ((m_R).^(theta - 1.0)) .- cost_level .* cost_elast .* (m_R.^(cost_elast-1.0))) ) .- 1 
    end 

    # for lower bound of root finding: use something close to zero 
    lower_bound = 1e-32 
    
    #### Step 2: Solve for m_R* using FOCs ####

    # upper_bound: can pick this as max (never want to spend beyond max subsidy (without constraints))
    upper_bound = ((theta .* epsilon_grid) ./ (cost_level .* cost_elast)).^(1.0 ./ (cost_elast - theta))

    # save results 
    optimal_m_R = zeros(length(epsilon_grid))

    # find root 
    @inbounds for row = 1:length(optimal_m_R)
        if epsilon_grid[row] < cost_level 
            optimal_m_R[row] = 0.0 
        else 
            optimal_m_R[row] = Roots.find_zero(m_R -> rent_seeking_FOC(m_R,z_tilde_grid[row], epsilon_grid[row]), (lower_bound,upper_bound[row]), Bisection())
        end 
    end 

    #### Step 3: Solve for remaining objects needed ####

    optimal_subsidy = epsilon_grid .* (optimal_m_R.^theta) .- cost_level .* (optimal_m_R.^cost_elast)
    optimal_revenue = z_tilde_grid .* ( (1.0 .+ optimal_subsidy ).^(param.elasticity_ratio .+ 1.0) )
    optimal_profits_C = (1-profit_tax) .* ((1.0 .- param.eta_tilde) .* optimal_revenue .- optimal_m_R)
    optimal_profits_NC = (1-profit_tax) .* (1.0 .- param.eta_tilde) .* z_tilde_grid
    TFPQ_model = (1.0 .+ optimal_subsidy) .* z_star_grid
    choose_C = (optimal_profits_C .- ((1-profit_tax) .* fixed_cost)) .>= optimal_profits_NC
    optimal_subsidy = ifelse.(choose_C,optimal_subsidy,0.0)
    optimal_m_R = ifelse.(choose_C, optimal_m_R, 0.0)
    optimal_revenue = ifelse.(choose_C, optimal_revenue, z_tilde_grid)
    optimal_profits = ifelse.(choose_C, optimal_profits_C .- ((1-profit_tax).* fixed_cost), optimal_profits_NC)

    #### Return final results ####
    return (
        z_star = z_star_grid, epsilon = epsilon_grid, z_tilde = z_tilde_grid,
        choose_C = choose_C, optimal_subsidy = optimal_subsidy, optimal_m_R = optimal_m_R,
        optimal_revenue = optimal_revenue, optimal_profits_C = optimal_profits_C, 
        optimal_profits_NC = optimal_profits_NC, optimal_profits = optimal_profits, 
        TFPQ_model = (1.0 .+ optimal_subsidy) .* z_star_grid, 
        revenue_output_C = z_tilde_grid .* ((1.0 .+ optimal_subsidy).^param.elasticity_ratio),
        revenue_output_NC = z_tilde_grid)
end 

# compute profits drawing epsilon grid
function compute_firm_problem_grid_industry(; z_star_grid, n_eps, param)

    # length 
    n_rows = length(z_star_grid)

    # create results vectors 
    expected_profits_C = zeros(n_rows)
    expected_profits_NC = zeros(n_rows)
    expected_m_R_C = zeros(n_rows)
    expected_subsidies_C = zeros(n_rows)
    expected_revenue_output_C = zeros(n_rows)
    expected_revenue_output_NC = zeros(n_rows)
    expected_subsidy_rate_C = zeros(n_rows)
    expected_revenue_C = zeros(n_rows)
    expected_expenditure_C = zeros(n_rows)
    subsidy_variation = zeros(n_rows*n_eps)
    epsilon_variation = zeros(n_rows*n_eps)
    epsilon_proba_variation = zeros(n_rows*n_eps)
    rent_seeking_share_variation = zeros(n_rows*n_eps)

    # for each row 
    @inbounds for row = 1:n_rows

        ### Substep 1: Draw epsilon conditional on z_star & the probability 
        Random.seed!(12345)
        normal_distrib = Normal(param.alpha_eps + param.beta_eps*log(z_star_grid[row]), sqrt(param.variance_eps_z))
        epsilon = rand(normal_distrib, n_eps)
        proba_epsilon = pdf.(normal_distrib, epsilon)

        # normalize probability to make sure it sums to 1  
        proba_epsilon = proba_epsilon ./ sum(proba_epsilon)

        if minimum(epsilon) <= 0.0
            println("Warning: Minimum epsilon is below zero: ", minimum(epsilon), ". Replacing by close to zero.")
            epsilon[epsilon .<= 1e-16] .= 1e-16
        end 

        ### Substep 2: Get optimal profits 

        results = compute_subsidy_industry(
            z_star_grid = ones(n_eps) .* z_star_grid[row], 
            epsilon_grid = epsilon, 
            param = param)

        # get expected profits by taking weighted mean across all possible epsilon 
        expected_profits_C[row] = sum(results.optimal_profits_C .* proba_epsilon)
        expected_profits_NC[row] = sum(results.optimal_profits_NC .* proba_epsilon)
  
        # collect further objects needed for later 
        expected_m_R_C[row] = sum(results.optimal_m_R .* proba_epsilon)
        expected_subsidies_C[row] = sum(results.optimal_revenue .* (results.optimal_subsidy ./ (1.0 .+ results.optimal_subsidy)) .* proba_epsilon)
        expected_revenue_output_C[row] = sum(results.revenue_output_C .* proba_epsilon)
        expected_revenue_output_NC[row] = sum(results.revenue_output_NC .* proba_epsilon)
        expected_subsidy_rate_C[row] = sum(results.optimal_subsidy .* proba_epsilon)
        expected_revenue_C[row] = sum(results.optimal_revenue .* proba_epsilon)
        expected_expenditure_C[row] = sum((results.optimal_revenue ./ (1 .+ results.optimal_subsidy)) .* proba_epsilon)

        # also save all subsidy variation 
        subsidy_variation[((row-1)*n_eps + 1):(row*n_eps)] = results.optimal_subsidy
        epsilon_variation[((row-1)*n_eps + 1):(row*n_eps)] = epsilon
        epsilon_proba_variation[((row-1)*n_eps + 1):(row*n_eps)] = proba_epsilon
        rent_seeking_share_variation[((row-1)*n_eps + 1):(row*n_eps)] = results.optimal_m_R ./ (results.optimal_m_R .+ (param.gamma_tilde .* results.optimal_revenue))
    end 

    expected_profits = param.proba_c .* expected_profits_C .+ (1.0 - param.proba_c) .* expected_profits_NC

    # adjust all variables by proba_c ? 
    expected_subsidies = param.proba_c .* expected_subsidies_C
    expected_revenue_output = param.proba_c .* expected_revenue_output_C .+ (1.0 - param.proba_c) .* expected_revenue_output_NC
    expected_revenue = param.proba_c .* expected_revenue_C .+ (1.0 - param.proba_c) .* expected_revenue_output_NC 
    expected_expenditure = param.proba_c .* expected_expenditure_C .+ (1.0 - param.proba_c) .* expected_revenue_output_NC 
    expected_m_R = param.proba_c .* expected_m_R_C

    return (
        expected_profits_C = expected_profits_C,
        expected_profits_NC = expected_profits_NC,
        expected_profits = expected_profits, 
        expected_m_R_C = expected_m_R_C,
        expected_m_R = expected_m_R, 
        expected_subsidies_C = expected_subsidies_C,
        expected_subsidies = expected_subsidies, 
        expected_revenue_output_C = expected_revenue_output_C,
        expected_revenue_output_NC = expected_revenue_output_NC,
        expected_revenue_output = expected_revenue_output, 
        expected_subsidy_rate_C = expected_subsidy_rate_C,
        expected_revenue_C = expected_revenue_C,
        expected_revenue_NC = expected_revenue_output_NC, 
        expected_revenue = expected_revenue,
        expected_expenditure = expected_expenditure, 
        expected_labor = (expected_revenue .* param.beta_tilde .* (1-vat_tax)),
        expected_capital = ((expected_revenue .* param.alpha_tilde .* (1-vat_tax)) ./ (rental_rate - 1.0)),
        subsidy_variation = subsidy_variation, 
        epsilon_variation = epsilon_variation, 
        epsilon_proba_variation = epsilon_proba_variation,
        rent_seeking_share_variation = rent_seeking_share_variation
        )
end

# Need function that solves across all industries & computes aggregates 
function solve_industries_baseline(; objects, n_eps = 500, io_elast = industry_io_elasticities)

    # create results object 
    results_firm_decisions = Vector{Union{Nothing, NamedTuple}}(undef, n_sectors)
    sectoral_expenditure = zeros(n_sectors) 
    sectoral_labor_demand = zeros(n_sectors) 
    sectoral_revenue = zeros(n_sectors) 
    sectoral_capital = zeros(n_sectors) 
    sectoral_rent_seeking = zeros(n_sectors) 
    sectoral_cit_revenue = zeros(n_sectors)
    sectoral_vat_revenue = zeros(n_sectors)
    sectoral_total_subsidies = zeros(n_sectors)
    sectoral_total_profits = zeros(n_sectors)

    # Compute profits grid for each sector (can parallelize this!)
    for sector in 1:n_sectors

        println("Sector ", sector, " of ", n_sectors)

        # Solve sector-level profits 
        results_firm_decisions[sector] = compute_firm_problem_grid_industry(
            z_star_grid = objects.z_star_grid[:,sector],
            n_eps = n_eps, 
            param = objects.param[sector])

        # Compute aggregate sector-level P_jY_j by aggregating p_ij y_ij (requires knowledge of distribution)
        sectoral_expenditure[sector] = sum( results_firm_decisions[sector].expected_expenditure .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry )

        # Aggregate sector-level labor demand 
        sectoral_labor_demand[sector] = sum( results_firm_decisions[sector].expected_labor .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry )

        # Also aggregate sector level revenue 
        sectoral_revenue[sector] = sum( results_firm_decisions[sector].expected_revenue .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry )

        # Aggregate other sector-level variables 
        sectoral_capital[sector] = sum( results_firm_decisions[sector].expected_capital .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry ) 
        sectoral_rent_seeking[sector] = sum( results_firm_decisions[sector].expected_m_R .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry )
        sectoral_cit_revenue[sector] = (profit_tax/(1-profit_tax)) * sum( results_firm_decisions[sector].expected_profits .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry )
        sectoral_vat_revenue[sector] = vat_tax * sum( ( (1-objects.param[sector].gamma_tilde) .* results_firm_decisions[sector].expected_revenue .- results_firm_decisions[sector].expected_m_R) .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry)
        sectoral_total_subsidies[sector] = sum( results_firm_decisions[sector].expected_subsidies .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry )
        sectoral_total_profits[sector] = sum( results_firm_decisions[sector].expected_profits .* objects.z_star_distrib[:,sector] .* objects.param[sector].share_industry ) 

    end 

    ## Construct (sectoral) aggregates 
    
    # aggregate labor demand 
    aggregate_labor_demand = sum(sectoral_labor_demand) 

    # get gamma_tildes 
    gamma_tildes = zeros(n_sectors)
    for sector in 1:n_sectors
        gamma_tildes[sector] = objects.param[sector].gamma_tilde 
    end 

    # construct Y_F (GDP)
    aggregate_final_use = sum(sectoral_expenditure) - sum( sum( io_elast .* gamma_tildes .* sectoral_revenue, dims = 1)' )
      
    ## normalize prices 

    # normalize initial relative prices 
    initial_sectoral_prices = ones(n_sectors) .* prod(network_nu .^ network_nu)

    # Solve for initial sectoral intermediate prices 
    initial_sectoral_interm_prices = (1 ./ prod( (io_elast[:,:]./initial_sectoral_prices).^io_elast[:,:], dims = 2))[:,1]

    ## Get further aggregates 

    # sectoral quantities 
    sectoral_output = sectoral_expenditure ./ initial_sectoral_prices

    # aggregate capital 
    aggregate_capital_demand = sum(sectoral_capital)

    # aggregate rent-seeking
    aggregate_rent_seeking = sum(sectoral_rent_seeking)

    # HH consumption under closed economy 
    total_HH_consumption_closed = aggregate_final_use - delta*aggregate_capital_demand - aggregate_rent_seeking # (enforce aggregate resource constraint -- only true if closed economy or BOP = 0) 

    # total tax revenues and subsidies 
    cit_revenue = sum( sectoral_cit_revenue ) 
    vat_revenue = sum( sectoral_vat_revenue )
    gross_tax_revenues = cit_revenue + vat_revenue
    total_subsidies = sum( sectoral_total_subsidies )
    net_govt_transfers = gross_tax_revenues - total_subsidies     

    # compute total HH income (not including capital)
    total_net_profits = sum( sectoral_total_profits )
    total_HH_income = net_govt_transfers + total_net_profits + aggregate_labor_demand # (wage is normalized to 1)
    total_HH_income_noprofits = net_govt_transfers + aggregate_labor_demand

    # report GDP net of rent-seeking 
    aggregate_final_use_net = aggregate_final_use - aggregate_rent_seeking

    ### return results ### 
    return(
        results_firm_decisions = results_firm_decisions, 
        sectoral_expenditure = sectoral_expenditure, 
        sectoral_labor_demand = sectoral_labor_demand,
        aggregate_labor_demand = aggregate_labor_demand,
        sectoral_revenue = sectoral_revenue, 
        aggregate_final_use = aggregate_final_use, 
        aggregate_capital_demand = aggregate_capital_demand, 
        sectoral_prices = initial_sectoral_prices,
        sectoral_interm_prices = initial_sectoral_interm_prices,
        sectoral_output = sectoral_output, 
        total_HH_income = total_HH_income, 
        total_HH_income_noprofits = total_HH_income_noprofits,
        total_HH_consumption_closed = total_HH_consumption_closed, 
        total_subsidies = total_subsidies,
        net_govt_transfers = net_govt_transfers, 
        gross_tax_revenues = gross_tax_revenues, 
        aggregate_rent_seeking = aggregate_rent_seeking, 
        sectoral_rent_seeking = sectoral_rent_seeking,
        sectoral_total_profits = sectoral_total_profits,
        total_profits = total_net_profits, 
        initial_sectoral_interm_prices = initial_sectoral_interm_prices,
        initial_sectoral_prices = initial_sectoral_prices,
        sectoral_subsidies = sectoral_total_subsidies,
        aggregate_final_use_net = aggregate_final_use_net
        )
end 
